rm(list=ls())
reference_system<-'+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs '

#Loading packages 

library(pacman)
p_load(sf,raster,dplyr)

#Path definition 

if(toupper(Sys.info()[['user']])=="MINERIA.RA"){
  
  # WS
  
  main_path<-"C:/Users/Mineria.ra/Dropbox/Research/PES Chiapas RCT/"

  
}else if(toupper(Sys.info()[['user']])=="SANTIAGO.FERNANDEZ"){
  
  # Santiago F UR
  
  main_path<-"C:/Users/santiago.fernandez/Dropbox/RA/PES Chiapas RCT/CreatedData"
  

}else if(toupper(Sys.info()[['user']])=="LENOVO"){
  
  # Santiago F Personal
  
  main_path<-"C:/Users/LENOVO/Dropbox/RA/PES Chiapas RCT/CreatedData"

}



setwd(main_path)
master<- read.csv("Lacandona_contracts.csv",sep=";")%>%
  mutate(Total_Property_folder_1=gsub("\\\\","/",Total_Property_folder_1))

reference_system<- '+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs'
months <- c("August2022","July2021","July2022","June2020","May2021" ,"May2022","Sept2022")

rasterDummy <- function(raster_path,shape_path){
  shape<-st_read(shape_path,quiet=T) ## Reading files
  raster<-brick(raster_path)
  shape<-st_transform(shape,crs=crs(raster)) #Get shape into raster crs
  addpoly <- as(shape, "Spatial") #Create polygon over raster
  plot(shape$geometry,add=T,border=color)
  new_raster<-raster(nrows=nrow(raster),ncols=ncol(raster),crs=crs(raster))
  new_raster<-setExtent(new_raster,extent(raster)) ##New raster with 1-0 from shape
  rasterized<-rasterize(addpoly,new_raster,1,getCover=TRUE)
  dummy_col<-values(rasterized)
  return(dummy_col)
}

#mapview::mapView(raster)+ mapview::mapView(shape)

for (i in 1:nrow(master)){
## Sample Raster to Rasterize
  raster_path <- paste0("Predictions/Final_Rasters/PropertyShapes/April2021/",
                        master$Polygon_reference[i],".tif")

  raster<-brick(raster_path)
  
  names(raster)[6:nlayers(raster)] <- months

  
  ## Property shape
   
  property_path <- paste0("PropertyContracts/PropertyShapes/",
                          master$Polygon_file[i])
  
  Contract_path <- paste0("Conafor_2021/Final_shapes/Contracted/",
                        master$Polygon_file[i])

  shape_path <- Contract_path

    if(master$PSA_application_file[i]!=""){
      
## CONAFOR file, not every shape has
      conafor_path <- paste0(
        "Conafor_2021/",master$PSA_application_folder[i],"/",
        master$PSA_application_file[i],".shp")
      }else{
        conafor_path <-   Contract_path
}
      ## Create Dummy by shape type
      color="blue"
      dummy_conafor <- rasterDummy(raster_path = raster_path,shape_path = conafor_path)
      color="red"
      dummy_contract <- rasterDummy(raster_path = raster_path,shape_path = Contract_path)
      color="green"
      dummy_property <- rasterDummy(raster_path = raster_path,shape_path = property_path)
      raster$Property <- dummy_property
      raster$Contract <- dummy_contract
      raster$Conafor <- dummy_conafor
      if(sum(values(raster$Contract))==0){
      print(i)
      }
        if (master$Contract_property_file[i]==""){
        raster$Contract <- 0
        }
        if (master$PSA_application_file[i]==""){
        raster$Conafor <- 0
        }

    writeRaster(raster,
            filename = paste0("pixel_id/Rasters/",
                              master$Polygon_reference[i],".tif"),
            format="GTiff",
            overwrite=T
            )

}
